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Abstract 

In the formalism of constrained mechanics, such as that which underhes the SHAKE and RATTLE meth- 
ods of molecular dynamics, we present an algorithm to convert any one-step integration method to a 
variational integrator of the same order. The one-step method is arbitrary, and the conversion can be 
automated, resulting in a powerful and flexible approach to the generation of novel variational integrators 
with arbitrary order. 



1 Introduction 

Consider a Lagrangian system defined by configuration space Q= {g} = M^, velocity phase space TQ = 
= M^^, and Lagrangian L: Q ^ M. Assume there is a holonomic constraint g{q) — 0, g: Q — > K'', 
suppose g has full rank, and let Q = g^^{0). The system evolves along curves q{t) e Q that are critical 
points of the action 

S= f L{q{t),v{t))dt (1.1) 

subject to the fixed endpoint constraints q{a), q{b) constant, and the first order constraint v{t) = q'{t). This 
variational principle is equivalent to the Euler-Lagrange equations 

where A(q, v) is found by solving the linear (Lagrange multiplier) problem (implicit sum on repeated indices) 



dv^dv^ " dq^ dv^dq^ dq^ ' , , 



—A" — 
dq^ dq^dqi 



V V 



These are the general Lagrangian systems. For example, they specialize to the Euler equations for the 
motion of a rigid body [TJ [12] by taking Q to be the 3x3 matrices {A} and to the the upper triangular 
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entries of AA^ — 1. And they specialize to the Kirchhoff approximations for the motion of an underwater 
vehicle [^I15j. The same explicitly constrained formalism is exploited in the molecular dynamics algorithms 
SHAKE and RATTLE [21 [S US] . 



The objective here is symplectic and momentum-preserving simulation of (1.2 1. Such simulations may 
be systematically generated by discretizing the defining variational principles, as in [lTTI 1131 IT7j . and partic- 
ularly |10j . which is specific to the constrained formalism. 

The theory of variational integrators is elaborated in O Uj, based on geometric discretizations of the 
velocity phase space TQ, i.e., based on 

1. certain assignments of curve segments to each {q,v) e TQ, where {q,v) is tangent to Q; and 



2. a discrete Lagrangian that approximates St for t = h; St is defined to be the classical action (1.1 1 



with a = 0, b = t subject to (1.2), and satisfies 



^ = L{q\v^). (1.4) 
at 

The discrete variational principle is a finite-dimensional constrained optimization problem, in which the 
objective function is a sum of the discrete Lagrangian on sequences in TQ. If the curve segments associated 
with the elements of the sequence agree to order r with the exact evolution of the Lagrangian system and 
the discrete Lagrangian agrees to order r with the exact classical action, then the variational integrator is 
order r accurate |14j . 

By definition, any one-step numerical integration method of order r gives order r accurate solutions 



to Equations (1.2) and (1.4). Any such method can be used to provide the curve segments and discrete 
Lagrangian that are required to construct a variational integrator as outlined above. In this article we derive 
the variational integrator from the corresponding discrete Euler-Lagrange equations in terms of a one-step 



integrator of (1.2) and (1.4) 



2 Basic Algorithm 

Given a Lagrangian L, a constraint g, and a one-step numerical integrator of order r, which wc call the 
standard layer, for the initial-value problem 

dq' , dv"- ^ dSt ^, . 

g(0) = q, i;(0) = V, 5((0) = 0, 
the aim is to generate a symplectic integrator of the same order. Let the standard layer be represented by 

{q,v)^Rt{q,v) = {R\{q,v),R\{q,v),R'l{q,v)). 

Assume that the standard layer exactly preserves the constraint; this restriction will be lifted later. By 
differentiating 5(g) = 0, the space of vectors tangent to Q = (7^^(0) is 

TQ={w={q,v):Dg{q)v^^]. 

Differentiating again. 



d_ 
de 



Dg{q + eSq){v + eSv) = v'^ g{q) 6q + Dg{q) 5v, (2.1) 

e=0 
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and so 

TTQ = {5w= {{q,v),{6q,Sv)) ■.Dg{q)Sq^O, v'^ g{q) Sq + Dg{q) Sv ^ O} . 

The notation T^TQ denotes the vector space {Sw : {w,Sw) £ TTQ}. TQ and TTQ are of course the first 
and second tangent bundle of the constraint. The notation for D^g{q) is problematic; it is a bilinear form 
with values in R*^. The quantity v'^D'^g{q), where v e M^, denotes the d x N matrix v'^d^g"' /dq^dq^ (q). 
We now construct the symplectic layer by defining the following quantities. 

1. The bias: a pair of numbers 

-1 < a" < 0, < a+ < 1, 
such that — a~ — 1. 

2. The time step: a number h > 0. 

3. The maps and the discrete Lagrangian L^: the required curve segments are associated with each 
element {q,v) by t i— > Rl{q,v) and the discrete Lagrangian is obtained from {q,v) R^{q,v). is 
not used. The ends of the segments provide the maps 

d-{q,v) = Rl^_{q,v), d+{q,v) = Rl^+{q,v), 

and the discrete Lagrangian 

Lh{q,v) = R^^^+{q,v) - Rf^-{q,v). 

4. The time step of the symplectic layer: Given wi = {qi,vi) € TQ, solve the following discrete Euler— 
Lagrange equations for W2 = ((?2,W2) € TQ: 

DLh{wi) 5wi + DLh{w2) 5w2 = 0, 

for all 5wi,5w2 satisfying 

Dd^{wi) Swi = 0, Dd+{w2) Sw2 = 0, 
Dd'l{wi) 5wi = Ddy^{w2) Sw2, 
iwi,Swi) e TTQ, {w2,Sw2) e TTQ. 

It is not necessary that the same standard layer provides both and 9^. For example, using any method 
to construct 9^, the adjoint [6] of the same method to construct 9^, and the bias a~ = = \, one 
evidently obtains a self-adjoint symplectic layer. Self-adjoint methods respect time reversal in the sense that 
a negative time step exactly reverses the discrete evolution. Also, self-adjoint methods are necessarily of 
even order: an odd-order self-adjoint method in fact has the next higher (even) order of accuracy because 
its odd order truncation errors must equal their negatives. 

The symplectic layer corresponds to the finite-dimensional discrete variational principle of finding the 
critical points of the discrete action 

Shiw, w) = Lh{w) + Lh{w) 
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subject to the constraints 

d^{w) — constant, d^{w) — constant, 

g{q,v)^0, g{q,v) = 0, 
as discussed in [4J. 

By the general theory, the maps 9^,9^ spht TTQ — ker Ddf^ ©kerD9^; i.e., {w,6w) sphts as 
Sw — 6w^ + Sw^ , 

6w^ e ker Dd^{w), 5w~ € ker Dd'^{w). 

(In the second hue, the presence of plus with minus is intentional and conforms to the notation of pT .) The 
discrete Lagrange one-form is defined by 

6~j^^{w) Sw = —DLh{w) 6w~ , 

and the general theory assures that the symplectic layer is a symplectic integrator with respect to a;^, = 

In Lagrangian systems, the Noether theorem shows that the presence of continuous symmetry is equivalent 
to the presence of conserved momenta. For example, translational [rotational] symmetry implies conservation 
of linear [angular] momentum; see [ll[12] for the basic theory, some of which the discussion here must assume. 
The discrete Noether theorem provides the same symmetry-momentum equivalence for discrete Lagrangian 
systems: Suppose that a symmetry group G acts on Q, such that 

1. 5 is invariant; 

2. ,dj^ intertwine the lift of the action to TQ; and 

3. £/i is invariant. 

Let be the Lie algebra of Q, and let £, E Q- Then the symplectic layer preserves the discrete momentum 
defined by 

where ^ € g and is the infinitesimal generator of ^ at w. For example, Q could be a matrix Lie group 
that acts on by matrix multiplication, the exact Lagrangian L invariant, the constraint g invariant, and 
the standard layer an explicit Runge-Kutta method. Then the standard layer intertwines the action on TQ, 
and Lfi is invariant, and the symplectic layer will preserve the corresponding discrete momenta. 

Neither the discrete symplectic form nor the discrete momentum equals in general the continuous coun- 
terpart. For example, if the system is rotationally invariant then the symplectic layer need not preserve 
some familiar mechanical angular momentum. The discrete Lagrangian system has its own version of angu- 
lar momentum, which is near to the continuous one, but through its dependence on 9^, 9^, is in general a 
complicated function of {q,v). 
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3 Constrained algorithm 

We require the following standard lemma, which will justify the use of a variety of (Lagrange) multipliers. 

Lemma 3.1. Suppose that A: E ^¥ and /i: E ^ G are linear. Let fj, be onto. Then A{e) — for all e such 
that ii{e) = if and only if there is a X: G ^ F such that A — A/i. 

Proof. Suppose A is zero on ker fi. Then A drops to A: E/ ker /i ^ F. Also, /i drops to ft : E/ ker /i — > G, and 
this is a linear isomorphism since /i is onto. Set A = Ap,~^. Then, letting w.E —^ E/ker/.t be the quotient 
map, A^(e) = Ap,~^iJ,{e) = ATr{e) = A{e). Conversely, if A = X/i and /Lt(e) = then A{e) = A(^(e)) 
A(0) = 0. □ 

We develop the equations used in the algorithm incrementally in stages; see Figures [T] and [2] In the figures, 
the dimension counts for the equations and variables are at right. The stages are equivalent representations 
of the same algorithm, starting from the fundamental description of the algorithm in Stage 0. Only the 
equations appearing in Stage 4 are implemented and solved in practice. In going from Stage to Stage 4, 
we increase the number of equations to be solved from 5(A^ — d) to 12N + 5d. 

Stage 0: The fundamental algorithm, given directly on the constrained phase space, as in J^. The funda- 
mental algorithm is defined on Q = g~^(0), regarded as a submanifold of Q. The constraints 

Dd-{wi) Swi = 0, Dd+{w2) Sw2 = 0, 

are enforced with multipliers (in TQ) A^ and A+, both having dimension dim Q = N — d. The constraint 

Dd'^{wi) Swi = Dd'^{w2) Sw2, 

is enforced with a multiplier fj,. Setting Swi and dw2 alternately to zero gives the two equations (SO. 2) 
and (SO.l), respectively. The connecting constraint d'^{wi) = d'^{w2) translates unchanged to (SO. 3). 

Stage 1: Enforce the restriction Swi G T^.TQ by introducing multipliers. The restriction {5q,5v) = Sw G 
T^TQ is 

Dg{q) Sq ^ 0, v^D^g{q) Sq + Dg{q) Sv ^ 0. 

For each equation in (SO.l) and (SO. 2), there are two corresponding multipliers vi, 1^2, which are row vectors 
of length d. 

Stage 2: Disambiguate A^,A^,/Lt. The multipliers A+, A^,/i are ambiguous as row vectors in up to any 
vector orthogonal to the constraint. Specifying them into TQ disambiguates them. Because Q = ^""'^(O) 
and the rows of Dg span the orthogonal complement to the tangent space of Q, the multipliers should have 
zero dot product with the rows of Dg, e.g., XDg(q)^ — for a multiplier X at q € Q. 

Stage 3: Lift the restriction that the standard layer preserves the constraint. A numerical integrator of order 
r used in the standard layer will in general preserve the constraint g only to accuracy order r. We posit a 
map i : X E'* ^ such that 

1. L{q, 0) = q; and 

2. if g{q) = 0, then Im DgL{q,0) is a complement of ker Dg; i.e., lmDei{q,0) ® ker Dg{q) = M.^ for all 
qeQ. 
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stage 0: Given wi = {qi,vi) 6 TQ solve, for W2 = {q2,V2) £ TQ, the equations 

(50.1) {OLh{wi) = X-Dd^{wi)+iJ.nd+{wi))5wi \/5wi eT-^^^TQ 2{N - d) 

(50.2) {nLhiw2) = \+Dd+{w2) - tJ.DdJ^{w2))5w2 Vit^a e T^^TQ 2{N - d) 

(50.3) d+{wi) = d-{w2) N-d 

5N -5d 

Lagrange multipliers A~ , A+, 3(A' — d) 

time advanced state {q2,V2) G TQ 2{N — d) 

5N -5d 

Stage 1: 

(51.1) DLh{wi) = X- Dd;; (wi) + ,iDd+ (wi) + u- [Dg{qi),0] + [v^ D^g{qi), Dg{qi)] 2N 

(51.2) DLh{w2) = X+Dd+(w2)-tiDd^{w2) + i^t[D9{q2),0]+iy+[v^D^g{q2),Dg{q2)] 2N 

(51.3) d+(wi) = d^{w2) N-d 

5N - d 

Lagrange multipliers A~, A+, 3(A' — d) 

Lagrange multipliers ,1^1 , ^2 4(i 

time advanced state (q2,V2) G TQ 2{N — d) 

5N-d 

Stage 2: 

(52.1) DLh{wi) = X-nd-iwi) + tJ.Dd+{wi) + u^[Dgiqi),0] + ly^ [vl D^'giqi), Dg(qi)] 2N 

(52.2) DLh{w2) = X+Dd+{W2) - I^Dd-{w2) + u+ [Dg{q2), O] + v+ [vlD^g{q2), Dg{q2)] 2N 

(52.3) a+ (wi) = a- (wi) N-d 

(52.4) X-Dg{d-{wi)y =0, iJ.ng{qY=0, X+ Dg{d+ {w2)Y = M 

bN + 2d 

Lagrange multipliers A~ , A+ , /i 3A'^ 
Lagrange multipliers v'^ ,1^^ , id 
time advanced state {q2,V2) G TQ 2{N — d) 

5N + 2d 

Stage 3: 

(53.1) DLh{wi) = X-Dd-(wi) + fiiDd+{wi) + iy-[Dg(qi),0] + [v^ D^g{qi), Dgiqi)] 2N 

(33.2) DLh{w2) = X+Df)+{W2) - fi2Dd-iw2) + 4[Dg{q2),0]+u+[v^D^g{q2),Dg{q2)] 2N 

(53.3) {X-r = DP{d-{wi)f{X-r, (A+)^ = DP(a+(«;2))^(A+r 2iV 

(53.4) /i^f = r>P(a+(«;i))^At, A2 = •DP(a;:(«)2))^M 2Ar 

(33.5) g = Pa+(w)i), 9^(«'2) = 2N 

(33.6) gr=P9^{«'i), 4=f'9h(^2) 2Ar 

(33.7) 9(132) = 0, Dg{q2)v2 =0 2d 

(33.8) A-Dg(gf)^ = 0, nDg{qf = 0, X+Dg{q+)'^=0 3d 

12N + 5d 

Lagrange multipliers A^ , A+, ^, A+, A~ , /ti, /i2 7N 
Lagrange multipliers , ly^ > '^1 i '^2 ^ 
variables q^ , q, 3^ 
variable 5+ d 
time advanced state {q2,V2) 6 TQ 2A' 

12N + 5(i 



Figure 1: Stages 0-3. 
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stage 4: Given wi = {qi,vi) e TQ solve, for W2 = (?2,t'2) 6 TQ, the equations 

(54.1) q-=Vd-{wi) N 

(54.2) q = Vd+{wi) N 

(S4.3a) A-Dq9-(iui) + i/f05(gi)+/iir)ga+(w)i) + i/2"«Jr-D^9(9i)--DgI'h(«)i) =0 AT 

(S4.3b) A-Dg(gf =0 d 

(S4.3c) (A-r = DP(9^(«;i))^(A')^ AT 

(S4.4a) /iiD„9+(t«i) + i/-Dg(qi) + A-D^a^(«;i)-D^Lft(tui)=0 AT 

(S4.4b) lJ.Dg(qY = d 

(S4.4c) iij = nV{d+{wi)Yti N 

{S4.5a) A+D,a+{«>2) + - fi2Dgd-{w2) + u+v'^ 0'^g{g2) - DgLh{w2) =0 AT 

(S4.5b) A+r>g(g+)^ = d 

(S4.5c) (X+)'^ = DV{d+{w2)yi\+r AT 

(S4.6a) D^Lhiwi) - u+ Dg{q2) - A+D„9+(i(;2) + [i2Ddh{w2) = AT 

(S4.6b) Dg{q2)v2 =0 d 

(S4.6c) /i2 = -0^(^,7 («'2))'^M AT 

(S4.7a) d-(w2) = i{q,e+) N 

(S4.7b) 9(92) = 0, d 

(S4.8) Q+=Pa+(«;2) JV 

12N + 5d 

Lagrange multipliers A~ , A+ , ^, A+ , A~ ,jli,jl2 7Af 
Lagrange multipliers , 1^2 > '^1 > '^2' 

variables q^ , q, q^ 3N 

variable 2d 

time advanced state (92, ^^2) G TQ 2N 



12N + 5d 

(S4.3a') (A- +//) + i/j-Go + u^vlD^giqi) - hDqL{qi,vi) = 
(S4.3b') A-GJ = 

(S4.4a') (/ia-A~ + ha+ ij) + v~Go - hvj Mq - hao = 
(S4.4b') tJ-Go = 

(S4.5a') (A+ - At) + 1/+G0 + u'vjD'^giqi) - hDqL{qi,vi) = 
(S4.5b') A+GJ = 

(S4.6a') tojMo - u^Go + hao - \+ + /j, = 
(S4.6b') G0U2 = 

(S4.7a') (92 - 9) - + ^a+i;2 = 

(S4.7b') Go{q2-q)=0 



Figure 2: Stage 4: the derivatives in stage 3 are split into partial derivatives with respect to q and v, and 
the equations are reorganized, the approximates are correspondingly the bottom primed equations. 
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Define a map P by 



(3.1) 



see the below figure. The map P follows the constant 0-fibers of l to where they intersect with Q. The maps 
^h'^h defined by the (constraint preserving) standard layer above. Letting the unconstrained standard 
layer define 9^ , 9^ , we redefine , 9^ by 



Also, we introduce the new variables 

g = 9^(wi) = 9,7 (wz), 
with variable and equation count each of N 

Q 




Usually an explicit projection P is not available and has to be computed iteratively. For example, t and 
P may be naturally defined by 

P{q) = q, 

q = q + Dg{q)^e, G (3.2) 

giq) = 0. 

As 9 varies, this particular i{q,9) moves g G Q away from Q orthogonally; the reverse, obtained from P, 
projects q orthogonally to Q. 

The connecting equation (S2.3) must be assumed to be full rank into Q; hence its (linearly independent) 
equation count is N — d. There are actually N equations when the image is considered into M.^ , as it must 
be for computations, but d of those are redundant because u'2 G Q by assumption and 9^, 9^ preserve the 
constraints. In the second of (S3. 5), explicitly writing the projection using new variable 9'^ resolves this 
problem because 9~^ robustly moves q away from the constraint Q, i.e., locally linearly in a nondegenerate 
way. So Equation (S3. 5) replaces (half of) Equation (S2.3), with equation count N rather than N — d, while 
the number of variables increases by d because that is the count for 6'+ . 

Stage 4- Split (S3.1) and (S3. 2) into partial derivatives with respect to q and v; rearrange terms and group 
equations. The q and v partial derivatives of Equation (S3.1) give Equations (S4.3a) and (S4.4a), respec- 
tively. Similarly, partial derivatives of (S3. 2) give (S4.5a) and (S4.6a). We group the equations so they can 
(eventually) approximated by linear equations with the same coefficient matrices, as will be seen. 
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4 Implementation 



We present here a strategy for the Stage 4 computation. 



4.1 The vector field and its derivatives. 

We now speciahze to Lagrangians of the form 

L = \mij{q)v'v^ + aj{q)v^ - V{q). 



(4.1) 



This is the most general quadratic Lagrangian with configuration-dependent coefficients. As is easily verified, 
Equation (1.3 1 becomes 



rriij A' 



dV 



dg^ 
dqi 



A' = 



Q2ga 



1 

2 



dq^dqi ' 
dmu ^ drriik 
dq^ dq^ 
_ daj_ 
dq^ ' 



dq^ 



(4.2) 



These are all linear equations for A and A with coefficient matrix of the form 

M{q)= [m,j{q)\. 



M{q) 
'Dg{q) 



-Dg{qY 




(4.3) 



The algorithm requires the derivatives of the maps 9^ and as well as the derivative of Lh ■ Automatic 
differentiation [5] can be used to compute these by computing the derivative of the one-step method in the 
standard layer that defines them. Alternatively, by Lemma 4.1 of jB], the derivative of the standard layer 
Rt that uses a Runge-Kutta method may be computed as the same Runge-Kutta method applied to the 
equations of first variation. In this case, it is only required to determine the derivative of the vector field. 

To compute the derivative of A^ with respect to q™ and v^' 



dA^ 




dg" 


Qqrn 


dq^ 


dq^ 


dA^ 


dXa 


dg- 






dq' 



dr. 

dq 



db, 
dq 



^3 y3 



dq^dq^ 



we differentiate the first equation in (4.2 1: 



dm., 
dq"" 



-A^ 



dq'^dq'' 



and for the constraints 

dg" dA' d^g" 



dq" dq"^ 
dg" dA' 

dq^ dv™- 



= 2 



dq^dq^dq^' 

Q2ga 



dq'-dq^ 



-A' 



These resulting equations are linear in the required derivatives with coefficient matrix (4.3 1. 
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4.2 Fixed-point iteration 

A possible approach for solving the Stage 4 (implicit) equations is by a fixed-point iteration. To find solutions 
to a generic equation f{x) = 0, split f = fo — A/ such that, for all b, an explicit solution to the equation 
fo{x) = 6 is available. We call fo an approximate. If A/ is sufficiently small, then in a suitable neighbourhood 
of the solution, the iteration Xi+i = fQ^(^Af{xi)) converges to a solution of f{x) = 0: 

foix) = Af{x) = /o(x) - fix) ^ fix) = 0. 
The iteration is, given an initial iterate xq, 

n+i = foixi) - fixi), solve faixi+i) = n+i, 
or equivalently, after substituting foixi) — ri, 



ra = /o(a;o), 
ri+i = ri - fixi), 
solve /o(xi+i) = r,;+i. 



(4.4) 



This approach is useful for the Stage 4 equations because they are nonlinear and a good choice for /o is gen- 
erally available. In this way, the Stage 4 computation may be organized into its equations and corresponding 
approximates. The required solution is obtained by iteratively evaluating the equations themselves and then 



solving for the approximates using (4.4 1. 
4.3 Stage 4 approximates 

Because the time step is small, the various configurations qi, q, q2, etc., are all close. Let Gq, Mq, and oq, be 
approximations to Dgiq), Miq), and a(q) respectively, obtained by evaluation at some such configuration; 
e.g., the configuration qi is a likely candidate. 



Equations (S4.1), (S4.2), and (S4.8) are of the form of Equations (3.2), which, for sufficiently small h, 
can be effectively approximated by Taylor expansion of giq) = at q = g: 

iq-q)- Gie = 0, 

- Go(q ~q)~ giq) - 0. 

Equations (S4.3c), (S4.4c), (S4.5c), and (S4.6c) involve the derivative of P in expressions such as 

= DViq)X^. (4.5) 



Differentiating equations (3.2 1 gives 



5q = Sq + Dgiq)-56 + D\e^g)iq)5q, 
^9 = Dgiq) Sq. 



(4.6) 



The matrix DFiq) is obtained by discarding SO after the inverse of (4.6), with Sg 0, i.e., 
DP(g) = [l 0] 



D^9-g)iq) Dgiq)- 
Dgiq) 



-1 


1 


0" 











10 



and Equation (4.5) becomes 







1 0' 




1 + D\e^g){q) Dgiqy- 


-1 


A^" 












Dg{q) 








If we define a variable z and put it in place of the zero in the matrix at left, then we can put a unit matrix 
in the (2, 2) slot of the first matrix on the right side, and invert. The result is the hnear equation 



+ D^{0-g){q) Dg{qy 
Dg{q) 

This replaces (S4.3c), (S4.4c), (S4.5c), and (S4.6c), and one can use the approximate 





A^' 




A^" 




z 








" 1 


/-it' 
Uq 




A^' 




A^" 


Go 







z 








(4.7) 



(4.8) 



which is computationally equivalent to 



1 












—Go 







z 








The remaining approximates are driven by the basic data of the variational principle: 
DqLh ~ hDqL, DyLh ~ hv^Mo + hao, 
d'^ {q,v) ~ q + ha'^v, df^ {q,v) ^ q + ha^v. 

In order, the approximates for Equations (4.3a,b) obtained from (4.9 1, A^ sa A 

DqdJ^{wi) fa 1, Dqd^ (wi) fa , 
Dqd+{wi) w 1, fiiDqd+{wi) w n, 
ViDg{qi) w j/f Go = 0, 
resulting in 

A" + ^i+v^Gq + i^^v^D'^giqi) - hDqL{qi,Vi) = 0, 



(4.9) 

and jl f^ fi are as follows: 



These are equations (S4.3a') and (S4.3b'). Similarly one obtains the approximates (S4.4a'b')-(S4.6a'b'), 
noting however that in (S4.5a') there are the further approximations 

v}vJD'^g{qi) sa i/2"wfD^5((ji), DqL{q2,V2) f^ DqL{qi,vi). 

Equations (S4.7a) and (S4.7b) are 

d^{w2)-Dg{qr9+ ^q, .9(92) = 0, 

which have to be solved for 92 and 9^ . Taylor expanding the second equation at q, and using g{q) — 0, gives 
the approximate 

q2-q-Dg{qir0+ + ha+V2 = 0, 

Go(g2 - g) = 0, 

which are Equations (S4.7a'b'). 
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4.4 Stage 4 solution 



The approximates (S4.3a'b') through (S4.5a'b') are solvable for the multipliers 
A", ^, A+, i^f , i^^, i^^, 



and the linear equations are all of the form (4.8). Indeed, one adds a~ times (S4.3b') and a+ times (S4.4b'), 
and then that together with (S4.4a') can be solved for a^X^ + a^fi and ■ Then, and similarly, (S4.3a'b') 
can be solved for X~ + n and z^f. Together these give X~ ,11, i/^ ,1^2 . Because n is then known, (S4.5a') with 
(S4.5b') minus (S4.4b') provide A+ and i^^ . In the same way (S4.6a'b') may be used to update z/^ and V2- 
Finally, (S4.7a'b') are solved for q2 and 

A , fii, fJ,2, X'^ 



occurring in (S4.3c), (S4.4c), (S4.5c), and (S4.6c'), may all be updated using appropriate versions of (4.7l 
and its approximates. The entire procedure can then be iterated until the variables q2 and V2 are at a 
predetermined accuracy. 

If should be noted that the multipliers that impose the constraints g, i.e., i^j^ , , z^j^ , 1^2 are not unique. 
For example, doubling g results in halving these multipliers. Such multipliers are nonphysical, and con- 
vergence of the iteration of Stage 4 should not be bound to the convergence of the multipliers themselves. 
Rather, the degree of convergence can be determined from products such as U2Dg[qi), which generally have 
the physical meaning of force of constraint; i.e., they are added directly in the equations to quantities with 
a physical interpretation. 
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